# scp savepoints/savepoint_2023-05-15/controls.rds UBELIX:/storage/homefs/jr18s506/projects/multilevel_wbe/savepoints/savepoint_2023-05-15/.
# scp savepoints/savepoint_2023-05-15/ww1.rds UBELIX:/storage/homefs/jr18s506/projects/multilevel_wbe/savepoints/savepoint_2023-05-15/.
# scp UBELIX:/storage/homefs/jr18s506/projects/multilevel_wbe/savepoints/savepoint_2023-05-15/ma5.3.1.rds savepoints/savepoint_2023-05-15/. 
# scp UBELIX:/storage/homefs/jr18s506/projects/multilevel_wbe/savepoints/savepoint_2023-05-15/ma5.3.2.rds savepoints/savepoint_2023-05-15/. 
# scp UBELIX:/storage/homefs/jr18s506/projects/multilevel_wbe/savepoints/savepoint_2023-05-15/ma5.3.3.rds savepoints/savepoint_2023-05-15/. 
# scp UBELIX:/storage/homefs/jr18s506/projects/multilevel_wbe/savepoints/savepoint_2023-05-15/ma5.4.1.rds savepoints/savepoint_2023-05-15/. 
if(!exists("controls")) controls = readRDS(fs::path("../savepoints/savepoint_2023-05-15/controls.rds"))
source("setup.R")
ww1 = readRDS(fs::path("../",controls$savepoint,"ww1.rds"))
shapes = readRDS(fs::path("../",controls$savepoint,"shapes.rds"))

We now model all ARAs together.

Switzerland

ww_all = ww1 %>%
  # log
  dplyr::mutate(logvl=log(vl)) %>%
  mutate(vl=if_else(vl==0, 1, vl)) %>%
  # create indexes for INLA
  dplyr::mutate(day1=day,
                ara1=as.numeric(as.factor(ara_n)),
                ara2=ara1) %>% 
  # group KLBS and KLZH (otherwise not identifiable)
  dplyr::mutate(lab2=if_else(lab=="KLBS","KLZH",lab),
                lab_n2=as.numeric(as.factor(lab2))) %>% 
  # group lab and method
  dplyr::mutate(lab_method=paste0(lab2,"_",method),
                lab_method_n=as.numeric(as.factor(lab_method)))

# correspondence table
corr_all = ww_all %>% 
  group_by(ara_n,ara_id,ara_name,kt,pop,lab,lab_n,lab2,lab_n2,lab_method,lab_method_n,ara1,ara2,NUTS2_name) %>% 
  count() %>% 
  ungroup() 
corr_all_ara = ww_all %>% 
  group_by(ara1,ara_name,ara_id,kt,NUTS2_name) %>% 
  count() %>% 
  ungroup() 

if(!controls$rerun_models) {
  mw_100_desc_table(ww_all) %>%
    dplyr::mutate(across(everything(),as.character)) %>%
    tidyr::gather() %>%
    flextable::flextable(cwidth=c(4,4))
}

key

value

Number of ARAs

118

Number of laboratories

9

Number of laboratory methods

11

Number of measurements

20535

Measurements below LOQ

235

Measurements below LOD

97

First

2022-02-07

Last

2023-05-14

Median viral concentration [gc/L]

1e+05 (range: 0 to 8e+06)

Median flow [m3/day]

1e+04 (range: 3e+02 to 9e+05)

Median viral load [gc/day/100,000]

4e+12 (range: 1 to 7e+14)

Model A5.3.1: unique national trend

We directly apply model A5.3 to all ARAs. We assume one temporal trend for Switzerland, and each ARA is free to deviate from it independently.

if(controls$rerun_models) {
  ma5.3.1 = INLA::inla(vl ~ 1 +
                         f(below_loq,model="iid") +
                         f(below_lod,model="iid") +
                         f(day,model="rw2", scale.model=TRUE, constr=TRUE,
                           hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
                         f(weekend,model="linear",mean.linear=0,prec.linear=.2) +
                         f(hol,model="linear",mean.linear=0,prec.linear=.2) +
                         f(ara1,model="iid") +
                         f(day1,model="rw1", scale.model=TRUE, constr=TRUE,
                           group=ara2, control.group=list(model="iid"),
                           hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))),
                       data = ww_all,
                       family = "gamma",
                       control.compute = list(waic=TRUE,config=TRUE),
                       control.predictor = list(compute=TRUE,link=1))
  saveRDS(ma5.3.1,file=paste0("../",controls$savepoint,"ma5.3.1.rds"))
} else {
  ma5.3.1 = readRDS(file=paste0("../",controls$savepoint,"ma5.3.1.rds"))
}
summary(ma5.3.1)
## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", " 
##    blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
## Time used:
##     Pre = 1.08, Running = 12194, Post = 12.8, Total = 12208 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) 13.925 7.710     -1.447   13.925     29.296 13.926   0
## weekend     -0.002 0.009     -0.020   -0.002      0.016 -0.002   0
## hol         -0.001 0.026     -0.053   -0.001      0.050 -0.001   0
## 
## Random effects:
##   Name     Model
##     below_loq IID model
##    below_lod IID model
##    day RW2 model
##    ara1 IID model
##    day1 RW1 model
## 
## Model hyperparameters:
##                                                 mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision parameter for the Gamma observations 3.448 0.048      3.351    3.450      3.538 3.456
## Precision for below_loq                        3.111 1.421      1.129    2.852      6.613 2.367
## Precision for below_lod                        0.012 0.006      0.004    0.011      0.029 0.009
## Precision for day                              0.011 0.003      0.008    0.011      0.018 0.010
## Precision for ara1                             1.064 0.253      0.702    1.017      1.683 0.916
## Precision for day1                             0.482 0.020      0.442    0.483      0.521 0.485
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1223613.19
## Effective number of parameters .................: 4231.92
## 
## Marginal log-Likelihood:  -732342.88 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
summary_exp_vl(ma5.3.1,pars="lab|method|hol|weekend")
##         exp(beta) 0.025quant 0.975quant
## weekend         1       0.98       1.02
## hol             1       0.95       1.05
ppp_vl_ara(ww_all,ma5.3.1)

avg_time_trend(ww_all,ma5.3.1)

mw_130_map_relative_vl(ma5.3.1,corr_all_ara,shapes)

mw_131_map_deviation_from_average(ma5.3.1,corr_all_ara,ww_all,shapes,12)

Model A5.3.2: effect of lab and method change

We add a covariate to measure the effect of the lab and of changes in methodology. To allow identifiability, we have to make sure that there are multiple ARAs per laboratory, so we group together KLBS (only 1 ARA) and KLZH (12 ARAs). The reference lab is ALTGR.

if(controls$rerun_models) {
  ma5.3.2 = INLA::inla(vl ~ 1 +
                         f(below_loq,model="iid") +
                         f(below_lod,model="iid") +
                         f(day,model="rw2", scale.model=TRUE, constr=TRUE,
                           hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
                         f(weekend,model="linear",mean.linear=0,prec.linear=.2) +
                         f(hol,model="linear",mean.linear=0,prec.linear=.2) +
                         f(ara1,model="iid") +
                         f(day1,model="rw1", scale.model=TRUE, constr=TRUE,
                           group=ara2, control.group=list(model="iid"),
                           hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
                         lab_method,
                       data = ww_all,
                       family = "gamma",
                       control.compute = list(waic=TRUE,config=TRUE),
                       control.predictor = list(compute=TRUE,link=1))
  saveRDS(ma5.3.2,file=paste0("../",controls$savepoint,"ma5.3.2.rds"))
} else {
  ma5.3.2 = readRDS(file=paste0("../",controls$savepoint,"ma5.3.2.rds"))
}
summary(ma5.3.2)
## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", " 
##    blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
## Time used:
##     Pre = 4.71, Running = 16149, Post = 4.68, Total = 16159 
## Fixed effects:
##                          mean       sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)            14.624 3015.142  -5898.535   14.624   5927.783 14.624   0
## lab_methodEAWAG_0       0.467    0.755     -1.012    0.467      1.947  0.467   0
## lab_methodEurofins_0   -3.076    0.574     -4.201   -3.076     -1.951 -3.076   0
## lab_methodEurofins_1   -1.522    0.577     -2.653   -1.522     -0.390 -1.522   0
## lab_methodKLZH_0        0.698    0.615     -0.509    0.698      1.904  0.698   0
## lab_methodLdU_0         0.019    0.654     -1.264    0.019      1.302  0.019   0
## lab_methodMicrosynth_0 -0.892    0.504     -1.880   -0.892      0.096 -0.892   0
## lab_methodMicrosynth_1 -1.262    1.694     -4.585   -1.262      2.061 -1.262   0
## lab_methodSCAV_NE_0    -1.125    0.980     -3.046   -1.125      0.797 -1.125   0
## lab_methodSUPSI_0      -0.713    0.809     -2.300   -0.713      0.874 -0.713   0
## weekend                -0.004    0.009     -0.022   -0.004      0.013 -0.004   0
## hol                    -0.007    0.025     -0.057   -0.007      0.042 -0.007   0
## 
## Random effects:
##   Name     Model
##     below_loq IID model
##    below_lod IID model
##    day RW2 model
##    ara1 IID model
##    day1 RW1 model
## 
## Model hyperparameters:
##                                                 mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision parameter for the Gamma observations 3.290 0.033      3.232    3.288      3.362 3.279
## Precision for below_loq                        1.579 0.553      0.601    1.519      2.690 1.448
## Precision for below_lod                        0.000 0.001      0.000    0.000      0.000 0.000
## Precision for day                              0.131 0.012      0.105    0.131      0.151 0.137
## Precision for ara1                             0.263 0.099      0.154    0.240      0.528 0.187
## Precision for day1                             0.512 0.017      0.477    0.513      0.544 0.515
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1223815.85
## Effective number of parameters .................: 4045.45
## 
## Marginal log-Likelihood:  -732451.58 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
summary_exp_vl(ma5.3.2,pars="lab|method|hol|weekend")
##                        exp(beta) 0.025quant 0.975quant
## lab_methodEAWAG_0           1.60       0.36       7.01
## lab_methodEurofins_0        0.05       0.01       0.14
## lab_methodEurofins_1        0.22       0.07       0.68
## lab_methodKLZH_0            2.01       0.60       6.71
## lab_methodLdU_0             1.02       0.28       3.67
## lab_methodMicrosynth_0      0.41       0.15       1.10
## lab_methodMicrosynth_1      0.28       0.01       7.85
## lab_methodSCAV_NE_0         0.32       0.05       2.22
## lab_methodSUPSI_0           0.49       0.10       2.40
## weekend                     1.00       0.98       1.01
## hol                         0.99       0.94       1.04
ppp_vl_ara(ww_all,ma5.3.2)

avg_time_trend(ww_all,ma5.3.2)

mw_130_map_relative_vl(ma5.3.2,corr_all_ara,shapes)

mw_131_map_deviation_from_average(ma5.3.2,corr_all_ara,ww_all,shapes,12)

Model A5.3.3: regional effects

Instead of just one temporal trend for Switzerland, we now allow independent temporal trends for each NUTS-2 region. ARAs- within each region are then allowed to deviate from the regional trend independently.

if(controls$rerun_models) {
  ma5.3.3 = INLA::inla(vl ~ 1 +
                         f(below_loq,model="iid") +
                         f(below_lod,model="iid") +
                         f(day,model="rw2", scale.model=TRUE, constr=TRUE,
                           group=NUTS2, control.group=list(model="iid"),
                           hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
                         f(weekend,model="linear",mean.linear=0,prec.linear=.2) +
                         f(hol,model="linear",mean.linear=0,prec.linear=.2) +
                         f(ara1,model="iid") +
                         f(day1,model="rw1", scale.model=TRUE, constr=TRUE,
                           group=ara2, control.group=list(model="iid"),
                           hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
                         lab_method,
                       data = ww_all,
                       family = "gamma",
                       control.compute = list(waic=TRUE,config=TRUE),
                       control.predictor = list(compute=TRUE,link=1))
  saveRDS(ma5.3.3,file=paste0("../",controls$savepoint,"ma5.3.3.rds"))
} else {
  ma5.3.3 = readRDS(file=paste0("../",controls$savepoint,"ma5.3.3.rds"))
}
summary(ma5.3.3)
## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", " 
##    blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
## Time used:
##     Pre = 3.02, Running = 10283, Post = 4.05, Total = 10290 
## Fixed effects:
##                          mean     sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)            14.766 16.104    -19.590   14.766     49.118 14.768   0
## lab_methodEAWAG_0       0.358  0.188     -0.010    0.358      0.730  0.357   0
## lab_methodEurofins_0   -3.035  0.181     -3.389   -3.035     -2.678 -3.036   0
## lab_methodEurofins_1   -1.838  0.196     -2.222   -1.839     -1.452 -1.839   0
## lab_methodKLZH_0        0.646  0.173      0.307    0.645      0.987  0.644   0
## lab_methodLdU_0        -0.148  0.175     -0.489   -0.149      0.197 -0.151   0
## lab_methodMicrosynth_0 -1.031  0.138     -1.300   -1.032     -0.760 -1.032   0
## lab_methodMicrosynth_1 -1.376  0.644     -2.639   -1.376     -0.113 -1.376   0
## lab_methodSCAV_NE_0    -1.078  0.265     -1.599   -1.078     -0.558 -1.077   0
## lab_methodSUPSI_0      -0.933  0.224     -1.372   -0.933     -0.492 -0.934   0
## weekend                -0.003  0.009     -0.020   -0.003      0.015 -0.003   0
## hol                     0.009  0.026     -0.042    0.009      0.059  0.009   0
## 
## Random effects:
##   Name     Model
##     below_loq IID model
##    below_lod IID model
##    day RW2 model
##    ara1 IID model
##    day1 RW1 model
## 
## Model hyperparameters:
##                                                 mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision parameter for the Gamma observations 3.301 0.043      3.223    3.298      3.391 3.290
## Precision for below_loq                        0.780 1.139      0.076    0.447      3.574 0.185
## Precision for below_lod                        0.004 0.004      0.000    0.002      0.015 0.001
## Precision for day                              0.005 0.000      0.004    0.005      0.006 0.005
## Precision for ara1                             9.554 1.243      7.164    9.537     12.048 9.630
## Precision for day1                             0.951 0.052      0.858    0.948      1.062 0.939
## 
## Watanabe-Akaike information criterion (WAIC) ...: 1223915.54
## Effective number of parameters .................: 3677.11
## 
## Marginal log-Likelihood:  -748599.68 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
summary_exp_vl(ma5.3.3,pars="lab|method|hol|weekend")
##                        exp(beta) 0.025quant 0.975quant
## lab_methodEAWAG_0           1.43       0.99       2.07
## lab_methodEurofins_0        0.05       0.03       0.07
## lab_methodEurofins_1        0.16       0.11       0.23
## lab_methodKLZH_0            1.91       1.36       2.68
## lab_methodLdU_0             0.86       0.61       1.22
## lab_methodMicrosynth_0      0.36       0.27       0.47
## lab_methodMicrosynth_1      0.25       0.07       0.89
## lab_methodSCAV_NE_0         0.34       0.20       0.57
## lab_methodSUPSI_0           0.39       0.25       0.61
## weekend                     1.00       0.98       1.01
## hol                         1.01       0.96       1.06
ppp_vl_ara(ww_all,ma5.3.3)

avg_time_trend_reg(ww_all,ma5.3.3)

mw_130_map_relative_vl(ma5.3.3,corr_all_ara,shapes)

mw_131_map_deviation_from_average(ma5.3.3,corr_all_ara,ww_all,shapes,12)

#’ ## Model A5.4.1: spatial structure #’ #’ #+ ma5.4.1a, fig.width=8, fig.height=12, R.options = list(width = 1000) if(controls\(rerun_models) { shapes_all = shapes\)ara_shp %>% dplyr::filter(ara_id %in% ww_all\(ara_id) %>% left_join(corr_all_ara,by = join_by(ara_id)) %>% dplyr::arrange(ara1) sf_use_s2(FALSE) graph_all = spdep::poly2nb(shapes_all) path_graph = paste0("../",controls\)savepoint,“W_all.adj”) nb2INLA(path_graph, graph_all) ma5.4.1 = INLA::inla(vl ~ 1 + f(below_loq,model=“iid”) + f(below_lod,model=“iid”) + f(day,model=“rw2”, scale.model=TRUE, constr=TRUE, group=NUTS2, control.group=list(model=“iid”), hyper=list(prec = list(prior = “pc.prec”, param = c(1, 0.01)))) + f(weekend,model=“linear”,mean.linear=0,prec.linear=.2) + f(hol,model=“linear”,mean.linear=0,prec.linear=.2) + f(ara1,model=“bym2”, graph=path_graph, scale.model = TRUE, constr = TRUE, hyper = list(theta1 = list(“PCprior”, c(1, 0.01)), theta2 = list(“PCprior”, c(0.5, 0.5)))) + f(day1,model=“rw1”, scale.model=TRUE, constr=TRUE, group=ara2, control.group=list(model=“iid”), hyper=list(prec = list(prior = “pc.prec”, param = c(1, 0.01)))) + lab_method, data = ww_all, family = “gamma”, control.compute = list(waic=TRUE,config=TRUE), control.predictor = list(compute=TRUE,link=1)) saveRDS(ma5.4.1,file=paste0(“../”,controls\(savepoint,"ma5.4.1.rds")) } else { ma5.4.1 = readRDS(file=paste0("../",controls\)savepoint,“ma5.4.1.rds”)) } summary(ma5.4.1) summary_exp_vl(ma5.4.1,pars=“lab|method|hol|weekend”) ppp_vl_ara(ww_all,ma5.4.1) #+ ma5.4.1b, fig.width=8, fig.height=6, R.options = list(width = 1000) mw_130_map_relative_vl(ma5.4.1,corr_all_ara,shapes) #+ ma5.4.1c, fig.width=8, fig.height=6, R.options = list(width = 1000) mw_131_map_deviation_from_average(ma5.4.1,corr_all_ara,ww_all,shapes,10)